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Abstract 

We perform a detailed study of the relaxation towards equilibrium in the Hamilto- 
nian Mean-Field (HMF) model, a prototype for long-range interactions in A^-particle 
dynamics. In particular, we point out the role played by the infinity of stationary 
states of the associated N ^ oo Vlasov dynamics. In this context, we derive a 
new general criterion for the stability of any spatially homogeneous distribution, 
and compare its analytical predictions with numerical simulations of the Hamilto- 
nian, finite N, dynamics. We then propose and verify numerically a scenario for 
the relaxation process, relying on the Vlasov equation. When starting from a non 
stationary or a Vlasov unstable stationary initial state, the system shows initially a 
rapid convergence towards a stable stationary state of the Vlasov equation via non 
stationary states: we characterize numerically this dynamical instability in the finite 
N system by introducing appropriate indicators. This first step of the evolution to- 
wards Boltzmann-Gibbs equilibrium is followed by a slow quasi-stationary process, 
that proceeds through different stable stationary states of the Vlasov equation. If 
the finite N system is initialized in a Vlasov stable homogenous state, it remains 
trapped in a quasi-stationary state for times that increase with the nontrivial power 
law N^''^ . Single particle momentum distributions in such a quasi-stationary regime 
do not have power-law tails, and hence cannot be fitted by the g-exponential distri- 
butions derived from Tsallis statistics. 
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1 Introduction 



Relaxation to tliermal Boltzmann-Gibbs equilibrium in A^-particle Hamiltonian sys- 
tems with long-range interactions has been recently the subject of an intense de- 
bate [1]. In some cases, the relaxation time has been shown to be extremely long and 
to increase with the number of particles. Hence, the study of these out-of-equilibrium 
conditions is of paramount importance for physical applications. The dynamics of 
systems with long-range interactions shows that some states are attained quickly, on 
time scales of 0(1), and that afterwards the system evolves very slowly, on time scales 
diverging with A^, towards Boltzmann-Gibbs equilibrium. We call states that evolve 
on time scales that diverge with N "quasi-stationary" . Some of them are character- 
ized by a particle distribution in the /x-space, /(r, p, t) = J^f S{r — rj(t), p — Pi{t)) 
(with (rj, pj) the position and conjugate momentum of the z-th particle and 6 the 
Dirac function), which remains close to a slowly varying smooth distribution for 
times that increase with A^. 

It should be remarked that most of the numerical evidences of this behavior are for 
ID and 2D systems^ . The theoretical explanation we propose in this paper, which is 
developed for mean-field models, extends to any dimension, as soon as the two body 
interaction is sufficiently smooth. Many recent studies of such quasi-stationary states 
are performed for the so-called Hamiltonian Mean Field (HMF) model (for a review 
see Ref. [3]). This model describes the motion of rotators under the action of an 
attractive or repulsive infinite range cosine interaction. In this paper we will consider 
the attractive case. The model then displays a second order phase transition, which 
is related to the development of a dynamical instability of the spatially homogeneous 
initial state with gaussian distribution of momenta, at a given value of the energy [4] . 

The analysis developed in this paper is based on a theorem due to Braun and 
Hepp [5,6], according to which the dynamics of a classical A^-particle system inter- 
acting via a two body and sufficiently regular long-range potential is well approx- 
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imated by the associated Vlasov dynamics of the density in /i-space. In the case 
in which the interaction force derives from a two-body smooth potential, Vlasov 
equation writes 

df 

^ + p-Vr/-Vrf/-Vp/ = . (1) 

The mean field macroscopic potential t/ is a functional of the probability distribu- 
tion function /(r, p, t), which makes the equation nonlinear in /. More precisely, 
Braun-Hepp's theorem states that, for a mean-field microscopic two-body smooth 
potential, the distance ^ between two initially close "weak" solutions of the Vlasov 
equation increases at most exponentially in time. The theorem applies also to singu- 
lar distributions, e.g. distributions that have a support on a set of dimension smaller 
than the one of the /i-space (for instance, a line in the two dimensional //-space of 
the HMF model). If we apply this result to a large particle approximation of 
a continuous distribution the error at t = is typically of order 1/\/N, thus for 
any "small" e and any "large enough" particle number A^, there is a time t up to 
which the dynamics of the original Hamiltonian and its Vlasov description coincide 
within an error bounded by e. The theorem implies that this time t increases at 
least as In A^. Extensions to wave particle dynamics of such a result have also been 
reported recently [7,8]. Since quasi-stationary states evolve on time scales that di- 
verge with A^, this result suggests that these states might gain their stability from 
being "close" to some stable stationary states of the Vlasov dynamics. 

Besides the conserved quantities of the Hamiltonian dynamics (energy, momentum, 
angular momentum...), the Vlasov description introduces additional integrals: the 
so-called Casimirs 

Cs[f] = Jsifiv,p,t))drdp, (2) 

for any function s. These conservation laws are responsible for the existence of an 
infinity of stable stationary states for the Vlasov dynamics. In this paper, we will 
argue that the existence of this infinity of stationary states is a possible explanation 
for the generic existence of quasi-stationary states in the finite A^-dynamics. 

This interpretation rises the following questions. May one predict the quasi-stationary 
states that emerge after a complex unstable Vlasov dynamics ? Among these sta- 
tionary states, are there some "statistically preferred" states ? What governs the 
relaxation towards equilibrium of the stationary states, and which is the scaling of 
the relaxation time with the number of particles? We will address in this paper the 
first and the third question, while the issues related to the second question will be 
only briefly mentioned. 

Before the slow relaxation phase sets in, a fast evolution takes place on a time scale 
that is independent of A^. This phenomenon is denoted as violent relaxation in the 
astrophysical context [9] , and is a consequence of the nonlinear complex dynamics of 
the Vlasov equation. After this very rapid stage, Vlasov dynamics produces thinner 
and thinner filament at ions of the density /, which lead to an apparent equilibrium 

^ The distance is measured in the Wasserstein metric, defined on the space of all measures. 
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described by a coarse grained density function /. A statistical mechanics interpreta- 
tion of this process has been proposed: for instance, for an initial two level density 
function, the equilibrium density is of Fermi-Dirac type [10]. Even if this may be 
already a partial answer to our second question above, we will not address in de- 
tail the very delicate point of the convergence of the Vlasov dynamics towards its 
statistical equilibrium in the present paper. 

Some authors [11] have advanced a challenging interpretation of the quasi-stationary 
states, suggesting that they should be true "equilibrium" states, obtained by max- 
imizing the Tsallis [12] entropy of the single particle distribution. We think that a 
landmark of such an interpretation would be the assessment of the existence of power 
law tails in the single particle momentum (or energy) distributions. Although we will 
not enter into a detailed fitting of such distributions using Tsallis g-exponentials, we 
will present a strong evidence of the absence of power law tails in the single particle 
momentum distribution at any stage of the time evolution, for the whole class of 
initial conditions we investigate in this paper, that are all homogeneous in space (see 
Section 5.3). For an initial condition in which the particles are concentrated in a 
point and momentum is uniform in an interval around zero, the authors of Ref. [11] 
have been able to fit a Tsallis g-exponential to the central part of the momentum 
distribution, after imposing an arbitrary cut-off to the tails. However, these authors 
do not exhibit any evidence of existence of power law tails, even for such a special 
initial state. 

Moreover, as mentioned above, our analysis applies also to this initial condition 
although it doesn't correspond to a stationary state of the Vlasov equation. At 
best, Tsallis statistics could describe the quasi-stationary states obtained from such 
a special initial state, but certainly not all of them, in particular those originated 
from the large class of homogeneous states that are studied in this paper. 

Moreover, the steadily progressing dynamical evolution observed for our class of 
initial conditions does not show any intermediate "statistically preferred" state. As 
stated above, the time evolution follows a sequence of stationary Vlasov states until, 
asymptotically, Boltzmann-Gibbs equilibrium is attained. 

It should be however mentioned that the special initial state studied in [11] is very 
interesting from a dynamical point of view, since it shows long-time correlations that 
are absent for the homogeneous initial states studied here (see also [13]). Similar 
initial states produce fractal structures in the /i-space for a self-gravitating sheet 
model [14]. 

As discussed above, Braun-Hepp's theorem suggests that the similarity between 
Hamiltonian A^-particle dynamics and Vlasov dynamics persists for times that in- 
crease as InA^. These times are linked with the fastest possible instability of the 
Vlasov dynamics. For stable solutions, the appropriate timescale is, however, the 
one associated with the fluctuations of the mean field. In agreement with this the- 
oretical remark, we present in Sec. 5.3 numerical results that indicate that the per- 
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sistence of quasi-stationary states is present up to times that are much longer than 
Braun-Hepp's In A^. This time scale increases as a power law in with a nontriv- 
ial exponent. Similar time scales have been found in gravitational systems. This is 
the case of Chandrasekhar's "coUisional" time scale, which is of order N/ IniV [15]. 
Although such time scale is similar to those we find in the HMF model, because it 
signals the final process of relaxation to Boltzmann-Gibbs equilibrium, its origin in 
our model is certainly different. 



In order to reformulate the timescale hierarchy sketched above, we can expect the 
following scenario: 



(1) An initial and fast evolution, well described for all initial conditions by the 
Vlasov dynamics, takes place on a timescale independent of the particle num- 
ber A^. 

(2) The system is then always trapped close to one of the numerous stable station- 
ary states of the Vlasov equation. This state may be the statistical equilibrium 
of the Vlasov equation (the most probable state with constraints given by the 
Vlasov invariants). 

(3) The system slowly evolves, on a much longer time scale, due to "collisions", 
or due to the fluctuations around this Vlasov stationary state. Consequently 
this time scale will be a function of A^. One can expect that this slow evolution 
takes place passing through different stable Vlasov stationary states. 

(4) Finally, the system reaches a particular Vlasov stable state, corresponding to 
the Boltzmann-Gibbs equilibrium state. 



We will try to give support to this scenario in the remaining of the paper. The plan is 
the following. We first introduce in Section 2 the Hamiltonian Mean Field model. We 
will then show in Section 3 that the Vlasov dynamics has an infinity of stationary 
states and we propose a new general stability criterion for any homogeneous dis- 
tribution (including non Boltzmann-Gibbs ones). We will compare these analytical 
results with numerical simulations of the finite- A^ Hamiltonian dynamics. The rapid 
convergence towards a stable stationary state of the Vlasov equation is described in 
Section 4. In the case of unstable stationary states, we show that the exponential 
destabilization may be investigated taking advantage of the existence of unstable 
modes of the Vlasov dynamics, in accordance with Braun-Hepp's theorem. The slow 
evolution towards equilibrium, passing through different stable Vlasov stationary 
states is described in Section 5. 
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2 The Hamiltonian Mean Field model 



We will consider the Hamiltonian Mean Field (HMF) model, whose Hamiltonian is 

j=l j,k=l 

where 6i G [— vr, 7r[ is the position (angle) of the z-th particle on a circle and pi 
the corresponding momentum. This system can be seen as representing particles 
moving on a unit circle interacting via an infinite range attractive cosine poten- 
tial, or as classical XY-rotators with infinite range ferromagnetic couplings [3]. The 
magnetization, defined as 

1 ^ 

M{t) = {M,,My) = -5:(cos%,sine,), (4) 

i=i 

or more precisely its modulus, M{t) = ||M(t)|| < 1, is the main observable that 
characterizes the dynamical and thermodynamical state of the system. 

Its introduction allows to write the canonical equations of motion as follows 



dt 
dpj 

'dt 



■Pj 



-M^sinej + My cos 6 j 



(5) 



Equilibrium statistical mechanics can be derived exactly both in the canonical and 
in the microcanonical ensembles [16,17,18]. In the ferromagnetic case, that we con- 
sider here, a minimal free energy (maximum entropy) state with a nonvanishing 
magnetization appears when lowering temperature (resp. energy per particle) below 
Tc = 1/2 (resp. Uc = 3/4). A discontinuity at Tc (JJc) in the second derivative of the 
free energy (entropy) with respect to magnetization signals a second order phase 
transition. This transition is between a low energy phase with particles forming a 
cluster (rotators pointing towards a preferred direction), and a high energy phase 
with particles evenly distributed on the circle (no preferred direction for rotators). 

This theoretical result, valid in the N ^ oo limit, is also confirmed by direct nu- 
merical simulations of the equations of motion (5), which moreover allows a careful 
analysis of finite corrections and give access to the study of non equilibrium 
features. It's in this context that, within the energy region U G [0.5, ?7c], a class 
of initial states has been found that displays an extremely slow relaxation towards 
Boltzmann-Gibbs equilibrium [11,16,19,20,21,22], with a relaxation time that in- 
creases with A^. Similar phenomena occur for other particle systems with long-range 
interaction (self-gravitating stars or point vortices [10]). Indeed, in the context of 
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two-dimensional fluid-dynamics [23,24] and plasma physics [24], the existence of an 
infinity of stationary states is known since a long time . 



These slow relaxation processes have recently attracted considerable attention, since 
the HMF model can be considered as a simple paradigmatic model of long-range 
interactions, without the two additional difficulties of gravitational dynamics: singu- 
larity at short range and particle evaporation. As briefly recalled in the Introduction, 
Latora et al. [11] have carefully analyzed an initial condition where all particles are 
located at the same position on the circle (giving initially M = 1) and momentum 
is uniformly distributed over a finite range, symmetrically around zero. The system 
shows a fast relaxation towards a small magnetization state which persists for an 
extremely long time, that increases with N . The authors compare the momentum 
distribution of such a quasi- stationary state with Tsallis distributions, obtaining 
some convincing fit of the central part of the distribution only after they impose a 
cut-off to momentum tails. Montemurro and Zanette [25], analyzing the same ini- 
tial condition, have even criticized the existence of a small magnetization plateau 
in time, by presenting some numerical evidence that magnetization first evolves to- 
wards a minimum, and then take off again towards the higher equilibrium value. 
We will avoid this controversial point by using the definition of quasi-stationary 
state given in the Introduction, i.e. we shall call quasi-stationary a state which still 
evolves, but on a time scale that diverges with A^. Hence, macroscopic properties are 
well defined over a sufficiently wide time span to allow local running time averages, 
even though the system slowly and continuously evolves towards equilibrium. 



At variance with most previous numerical experiments [11,25], we choose here an 
initial state where the particles are uniformly distributed on the circle (hence M is 
initially close to zero, M = 0(1 /\/N)) and momentum has a uniform distribution 
centered around zero, as above. We make this choice for two reasons: i) this state 
is a stationary state of the Vlasov equation, that describes the HMF model in the 
N oo limit, ii) the M = state plays a relevant role also in the previous numerical 
experiments [11,25] and we then thought that it is better to start directly from it. 
Finally, we think that the analysis of the M = initial state can also clarify several 
aspects of the phenomenology of the M = 1 initial state used by the authors above. 



As explained in the Introduction, we believe that the Vlasov description is a useful 
and appropriate tool to understand the slow relaxation process. In the following 
Section, we introduce and discuss the Vlasov equation for the HMF model. 
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3 The Vlasov dynamics of the HMF Model 



3. 1 Introduction 



In the continuum limit, that is keeping the volume (here the interval [— vr, 7r[) and the 
energy per particle fixed as the number of particles — > oo, the dynamics governed 
by Eqs. (5) is described by the Vlasov equation. The state of the finite system 
can be described by a single particle time- dependent density function 

1 ^ 

h{o.P.t) = -Y.He-e,{t),p-p,{t)) (6) 

where 5 is the Dirac function. When is large, it is natural to approximate the dis- 
crete density fd by a continuous one / t). Using this density, also called //-space 
distribution, it is possible to rewrite the two components of the magnetization M 
given by Eq. (4). They read 



MM= J fi^^P^t) cos^dMp , (7) 
My[f]= J f {9, p,t) sine dedp . (8) 

Within this approximation, one can write the potential that affects all the particles 

as 

(^) [/] = -M,[/] cos ^-M,[/] sine . (9) 
This potential enters the expression of the Vlasov equation 

which governs the spatiotemporal evolution of the density /. In the remaining of 
the paper, we will omit the over-bar on and My for the sake of simplicity. 



It is important to notice that the discrete distribution fd, a sum of Dirac peaks, 
contains exactly the true dynamics of the system and is also a solution of the Vlasov 
equation (10). Then, introducing a suitably defined distance on the space of proba- 
bility measures on [— vr, 7r[x]R, it is possible to show [5,6,7] that the distance between 
two solutions /i and /2 of Eq. (10) grows at most exponentially in time 

rf(/i(t),/2(t))<rf(/i(0),/2(0))e^^ , (11) 

with a growth rate 7 which is independent of the initial conditions. This result, which 
is the essence of Braun-Hepp's theorem, heavily relies on the mean field character 
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of the underlying Hamiltonian dynamics and on the genericity of exponential insta- 
bility of trajectories. Choosing then fi{9,p,t) = fd{9,p,t) and taking for f2{0,p,t) 
a continuous approximation of fd, one immediately obtains that d (/i(0), /2(0)) — >• 
as N grows. As a consequence, the previous result implies that the discrete particle 
dynamics converges to Vlasov dynamics when iV — *• oo, uniformly over all fixed time 
intervals [0, T]. However, for all fixed there is a typical time r ~ 7"^ over which 
the two dynamics diverge. 

From this analysis, it is therefore natural to expect that particle and Vlasov dynam- 
ics coincide during a time that diverges as InA^, if (i(/d(0), /(O)) ~ and 7 is 
independent of A^. However, for initial conditions corresponding to stable stationary 
solutions of the Vlasov equation, this time may be much longer, actually of order 
A^ [26]. 

All this explains why the properties of the Vlasov equation are of particular interest 
for the study of the particle dynamics. In the next subsection, we will study the 
Vlasov dynamics of the HMF model and its stationary states, with the aim of getting 
useful insights on particle dynamics. 

3.2 The Vlasov dynamics and its stationary states 

The Vlasov equation inherits from the particle dynamics the conservation of the 
energy Hy 

Hv[f] = l^f{0,P,t)dedp + ^--^^ (12) 
and of the total momentum 

P[f] = Jpf{9,p,t)dedp . (13) 

However, the Casimir's functionals (2), defined for any continuous function s, yield 
an infinity of additional conserved quantities, linked with the labeling symmetry 
when following a fluid particle in the /i-space. These new conserved quantities play 
of course a major role in the dynamics of the Vlasov equation and thus of the particle 
dynamics. 

For a non stationary initial distribution, the dynamics of the Vlasov equation is 
known to give rise to a very complex nonlinear evolution, characterized by stretch- 
ing and folding of the initial distribution, the details of this evolution being usu- 
ally unpredictable. However, one may predict the final evolution using statistical 
mechanics arguments, in the spirit of the statistical mechanics of two dimensional 
conservative flows [27,28] or of the Vlasov-Poisson equation [29,30]. One ends up 
with the most probable coarse-grained distribution /, which takes into account the 
dynamical invariants. Unfortunately, the dynamical mixing of the distribution is 
likely to be incomplete: the reason lies in the existence of infinitely many stable sta- 
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tionary solutions of the Vlasov equation, in the neighborhood of which the system 
may be trapped. Therefore, Vlasov dynamics quickly converges (weak convergence) 
to a stable stationary state which should be studied in detail. 

Equations (7) and (8) show that the magnetization M is constant for a stationary 
solution f{6,p). This implies that the potential V is constant. The equation for the 
stationary states of the Vlasov equation may thus be considered as a linear first 
order partial differential equation. Solutions are then given by densities / that are 
constant on the characteristics of the equation, corresponding to the level sets of the 
individual particle energy 

2 2 

eie,p) = ^ + Vie) = ^-M^cose-MySme , (14) 

which corresponds to the energy of a pendulum. It is important to observe that 
the total energy Hy given by Eq. (12), is different from the sum of the individual 
energies (14). 

Smooth stationary solutions of the Vlasov equation are thus given by 

f{e,p) = <^{e{9,p)) , (15) 

where $ is any real function. Moreover, the values of M^, My and the function $ 
must be self-consistent. The Boltzmann-Gibbs equilibrium density 

f,,i9,p)=Aexpi-Peie,p)) (16) 

is a particular case, although very important. One may also prove that stationary 
states of the Vlasov equation in a moving frame with constant velocity v are given 
by /(^,p) = $(e {9,p)+vp). 

Let us note that the function $ may be multi- valued for individual energies e greater 
than the energy of the separatrix (one branch for particles with positive momentum 
and the other for negative momenta). We will assume that this does not happen 
in the remaining of the paper, for the sake of simplicity (a generalization would be 
straightforward) . 

3. 3 Stability of stationary states of the Vlasov equation 

As discussed above, the stationary states of the Vlasov equation are not true sta- 
tionary states of the particle dynamics. If they are stable, however, they may explain 
the long lifetime of quasi-stationary states in the particle dynamics. Linear stability 
results for the stationary states of the Vlasov equation have been already reported 
in the case of spatially homogeneous distributions for both Gaussian [4] and water 
bag [16] momentum distributions. In this Section, we will show that it is possible to 



10 



derive stability results for arbitrary spatially homogeneous stationary states, using 
a method developed in the context of two dimensional fluid dynamics and in plasma 
physics [24], based on original ideas introduced by Arnold [23]. 

The authors of Ref. [24] actually study nonlinear stability, a concept that we would 
like to briefly distinguish from other stability concepts. For a generic dynamical 
system, any extremum /o of a conserved quantity F [/] is a stationary point of the 
dynamics. It is said to be formally stable if the second variations A2-F [5/i, 5/2] 
of F is positive definite (/o is then a minimum) or negative definite (/o is then a 
maximum). In the case of the linearized dynamics around a formally stable point 
/o, as the second variations of F at /o are conserved, a small perturbation of /o 
remains bounded in the norm provided by the second variations: this state is linearly 
stable. Since this implies that the spectrum of the linearized dynamics does not 
have any negative value, the system is also spectrally stable. It is however not true 
in general that spectral stability implies linear stability, and that linear stability 
implies formal stability. Finally, nonlinear stability corresponds to the case where 
a small perturbation, evolving according to the real dynamics, remains bounded in 
some norm. It can be shown that nonlinear stability implies spectral stability, the 
converse being wrong in general, whereas formal stability implies nonlinear stability 
only in finite dimensional systems. 

In this Section we prove that any stationary state of the Vlasov equation, defined by 
Eq. (15) with $ strictly decreasing, corresponds to a critical point of some invariant 
functional. Computing the second variations of this functional, we can therefore 
exhibit a necessary and sufficient condition of formal stability for such a stationary 
state. 

Let us consider the maximum of the functional 



where Hy is the energy (12), Cs is a Casimir functional (2) corresponding to a 
strictly concave function s and /5 is positive. Performing the first variations of this 
functional, we obtain the equation 



s' ifo) = (3 [ — - cos {6 - a) fo{a,p,t)dadp] + jj, = /3e {6, p) + fi , (18) 



with e given by (14) and \l/ the inverse function of s', the derivative of s. The condi- 
tion that s is strictly concave is equivalent to the fact that \1/ is strictly decreasing. 



F [/] = Cs [/] - PHv [f]-f^l f{e, p, t) dMp 



(17) 




which defines the critical points /q. This yields 



/o(^,p) = vl>(/5e(^,p)+/i) , 
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The computation of the second variation of F gives: 

A2F[6f,6f] = J s"{fo{e,p))[6f{e,p)]'dedp + (3{{M,[6f]f + {My[6f]f) , (19) 

where M^. and My are given by Eqs. (7) and (8). As s is strictly concave, s" is 
negative. The first term of A2F [6f, 6f] is thus clearly negative whereas the second 
one is positive. We will consider now only homogeneous states, corresponding to 
M, = My = and /o (O^p) = fo (p) = ^ W/2 + /i). 

Let us introduce the Fourier series of the perturbation 

5f {6,p,t) = '^Cn{p,t) cosnO + Sn{p,t) smnO . (20) 

n 

From Eq. (19), after integration on the spatial variable 6*, one obtains the second 
variations 

A2F[5/,5/]= /s"(/o(p))E(4(p) + 4(p))dp+2G(ci(p)) + 2G(.i(p)) , (21) 

n>l 

where we have introduced 

G (c (p)) = j s" {fo (p)) c' (p) dp + ^ (/ c (p) dp)' . (22) 

The terms involving c„ and s„, n > 1, in Eq. (21), are clearly negative definite, since 
s" is strictly negative. Consequently, the second variations of F are negative definite 
if and only if G is negative definite. 

The sign of the function G can be studied using the Cauchy-Schwartz inequality: 



e,rtdH^f/^n^^*| (23) 



v'-»" (A (p)) 



<(/."(/oW)c^(p)dp)(/^^). (24) 
This inequality leads therefore to 

G (c (p)) < / s" (/o (p)) (p) dp f 1 + ^ f / ,, dp] ] . (25) 



s"ifoip)) 

Recalling that s" is strictly negative, we conclude that if the quantity 

is positive, the function G is negative. Conversely, when it is negative, considering 
the particular function c (p) = —1/s" (/o (p)) demonstrates that G may be positive. 
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Differentiating witli respect to tlie variable p the equality 

^'(/o(p))=/5y + /i, (27) 
obtained from Eq. (18) in the case of homogeneous states, yields: 

s"{fo{p)) P ' ^ ' 

As /o(0) = 0, the quantity (26) can be written as 

I[fo] = l + l r^dp . (29) 

^ J — oo p 

This leads to the following equivalence: 

/o is formally stable <^==^ J[/o] > . (30) 

This condition will of course be an excellent criterion to predict the stability of 
several initial conditions, and therefore the lifetime of the corresponding quasi- 
stationary states. This is what we will consider now. 

Let us note that Inagaki and Konishi [31] have found a dispersion relation for the 
linearized Vlasov equation which leads to the above stability criterion. Hence, in 
this particular case, linear stability and formal stability criteria coincide. 



3.4 Applications of the nonlinear stability criteria 



3.4.0.1 Waterbag distribution One of the most widely used initial condition 
in previous numerical studies of the HMF model is the so called waterbag distribu- 
tion 

{0 if Ipl > p 

(31) 
1 / (2p) if — p < p < p . 

If Mj. = My = (homogeneous state), the relation between energy Hy = U and p 
is p = V6f/ — 3. Computing the first derivative of /„b 

fUip) = ^[KP + P)-S{p-p)] , (32) 



one obtains the following expression 



I[U] = l-~ , (33) 
z p 
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exhibiting that the critical width of the distribution, above which the waterbag is 
formally unstable, is pc = This corresponds to the critical energy density 

U: = ^ , (34) 

as reported earlier [16]. 



3.4.0.2 Gaussian distribution As a second example, let us consider a Gaus- 
sian distribution, 

/.(P) = /J e-^^^/^ , (35) 

characteristic of an equilibrium canonical distribution. In this case, the quantity (29) 
can be easily reduced to 

P 
2 

emphasizing that the critical inverse temperature is /?* = 2 and consequently the 
critical energy density U* = 3/4. This coincides with the critical point of the second 
order phase transition Uc- 



/[/g] = l-^ , (36) 



3.4.0.3 Mixed distribution Finally let us consider a more general distribution, 
namely a mixed distribution between /^b and /g, defined as 

f^(p) = (l-a)U{p) + af,{p) . (37) 

Thanks to the linearity of the quantity (29) with respect to the distribution, the crit- 
ical energy density for this mixed distribution fa is obtained as a linear combination 
of both previous results: 



t^:(a) = ^(l-«) + ^« = ^ + ^ • (38) 

Such a result allows to define a phase diagram of the dynamical critical energy that 
we will be able to confirm using numerical simulations. We aim in the following 
at showing how these considerations about the Vlasov equation and its stationary 
states are useful to understand the dynamics of the discrete particle system. We 
begin in Section 4 by studying the short time evolution, and turn in Section 5 to 
the intermediate and long time behaviour. 



14 



4 Short time behavior 



4-1 The numerical setup 



We have numerically integrated the canonical equations of motion (5), by using 
symplectic 4th-order integrators, the McLachlan-Atela's [32] or Yoshida's [33] algo- 
rithms. The timestep dt = 0.1 was chosen to obtain an energy conservation with a 
relative accuracy \ AE/E\ better than 10"''. We will consider initial conditions with 
uniform distribution with respect to 6 as explained in Section 3.3. The magneti- 
zation being consequently zero, these states are stationary solutions of the Vlasov 
equation (10). 

However, the numerical calculations correspond to simulations with a finite number 
of degrees of freedom, and these initial states are not anymore stationary a priori. 
The spatial coordinates 6j are randomly chosen in the interval [— 7r,7r[, which leads 
to a magnetization of order 1 / ^/N initially. The momenta pj are also randomly 
chosen from the given distribution f{p) satisfying the conditions for the energy 

±ytl = u-- , (39) 
N^-, 2 2' ^ ^ 

whereas the conserved total momentum is set to zero 

Ep.— O . (40) 

We will present how the numerical results allow us to detect the critical energy 
for dynamical instability, using first the most widely used initial conditions, the 
waterbag distribution, and then the mixed ones. 



4-2 Waterbag initial distributions 



4-2.1 The first peak of the magnetization 

Figure 1 presents the initial temporal evolution of the magnetization M(t) for differ- 
ent values of the energy density, but with the same number of particles, for a system 
initialized with a waterbag distribution in momentum and a homogeneous one for 
the angles(M is close to zero since initially M = 0{1/ \/N)). Averages over a set of 
initial conditions (samples) are taken. These results (already partially reported else- 
where [34]) show that the initial time evolution of the magnetization, starting from 
such a non equilibrium initial condition, is quite different depending on whether U is 
larger or smaller than U* ~ 0.583. If f/ > f/*, but still below the second order phase 
transition energy Uc, the magnetization remains close to the M = initial value 
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and does not show any tendency towards the nonvanishing equihbrium value (the 
long time relaxation to equilibrium will be discussed in Section 5.3). For U < U*, 
instead, the magnetization shows a fast relaxation to a non vanishing value which is 
close to equilibrium. Relaxation proceeds through repeated oscillations that damp 
after a relatively short time. In order to characterize quantitatively this behavior, let 
us focus on the first peak of M(t), by studying its height and its time of occurrence 
as a function of the energy density U, as presented in Fig. 2 for increasing particle 
numbers. 




Fig. 1. Temporal evolution of M{t) for different values of the energy density U, when the 
number of particles is = 10'^. The values of U are from 0.50 to 0.78 with 0.04 step size. 
The curves correspond to averages over 100 samples. 




35 
30 
25 
20 
15 

10 ; 

5 









(b) ■ 






< ^=10" 




' N=W^ ° , 























0.5 



0.55 



0.65 

u 



0.7 



0.75 



0.8 



Fig. 2. The first peak height (a) and the first peak time (b) vs. energy density U . The 
curves correspond to different particle numbers N = 10^,10^,10^ and 10^ from top to 
bottom in [/ > [/* for panel (a) and from bottom to top in ?7 < [/* for panel (b), keeping 
in mind that every curve is obtained averaging over 100 (resp. 20) samples for N < 10^ 
(resp. N = 10^). The vertical line represents the theoretical critical energy density U*, 
given by Eq. (34). 

Figure 2(a) emphasizes that the first peak height vanishes in the energy region 
above U* as N increases, in agreement with the theoretical predictions of the Vlasov 
equation, which imply that the M = state becomes stable in this energy range in 
the continuum (A^ — > oo) limit. The initial fluctuations of the magnetization, which 
are of order 1/\/N, do not grow when U > U*, because the waterbag state is Vlasov 
stable. On the contrary, these fluctuations rapidly grow for U < U*, leading in a 
short time to a nonvanishing magnetization state which is close to equilibrium. 
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Similar indications come from the first peak time, shown in Fig. 2(b). For U > U* 
this quantity clearly shows a convergence to an asymptotic value as increases, 
but this value sharply increases as one approaches U* from above, signaling the 
instability of the M = initial state. The behavior below U* is less clear since the 
first peak time does not yet show the saturation expected on the basis of the Vlasov 
equation for the values of considered here. This might perhaps be due to the 
limited validity in time of the Vlasov description when starting from an initially 
unstable state. We will however show in the next Section that at least the early 
exponential instability is well reproduced. 

Both quantities, the first peak height and the first peak time, are therefore useful 
indicators of the presence of a dynamical critical energy U* , as predicted theoretically 
in the continuum limit. To summarize, we could say that on a 0(1) time scale, and 
when starting from such a waterbag initial state, one would observe the instability 
of the M = state at U* instead of Uc, the phase transition energy. 

Before extending these findings to other types of initial distributions in Section 4.3, 
we will first show in the next Section that the predictions of the Vlasov equation 
are sufficiently precise to give the growth rate of the instability. 



4-2.2 Theoretical estimate of the growth rate 



As shown in Fig. 2(b), a typical time scale of the system (the one of the first peak 
in magnetization) shows a tendency to diverge at U*. We show in this Section 
that Vlasov equation not only predicts this divergence, but is also quantitative in 
determining the time scale of the initial exponential instability of the magnetization 
for U <U*. 

Let us therefore define the exponential growth rate A of an initial perturbation 
6f{6,p,0) of a Vlasov stationary unstable state as follows 



\\6f{9,p,t)\\-exp{Xt)\\6f{9,p,0)\\ . (41) 

The perturbation we consider is around a homogeneous state, whose density func- 
tion / depends only on the variable p. However, the perturbation 6f depends on all 
variables {9,p, t). 



The linearized equation that governs the time evolution of the perturbation 6f, for 
vanishing magnetization, is 

= ° • ^^^^ 

Using the development in Fourier series of the perturbation 6f given in Eq. (20), we 
obtain the following equations for each Fourier component 



17 



dCn 




dt ~ 


-npsn 


dSn 




dt 




dci 




'dt^ 


-psi - 


dsi 




'dt^ 


pci + 



Idl 

4 dp 
Idf 



Vn > 1 (43) 
Vn > 1 (44) 
du si{u,t) (45) 



-J- J du c,{u,t) . (46) 

One can easily show that the Fourier components with n > 1 cannot be unstable. 
Indeed, by introducing the quantity Vn = Cn + iSn, it is straightforward to show that 
its time derivative is v'^ = ipnvn- Hence, the generalized eigenvalues are all pure 
imaginary numbers, inp (with p G M), and the corresponding eigenvectors are Dirac 
delta functions. 

Let us then concentrate on the n = 1 components for waterbag initial conditions. 
Using expressions (32) for the derivative with respect to p of the distribution (31), 
Eqs. (45) and (46) become 



dc 1 

-psi - — [6{p + p) - 6{p - p)] I du si{u, t) (47) 



dt 8p 

^= pci + ^[6{p + p) -6{p-p)] J duci{u,t) . (48) 

To solve this infinite dimensional eigenvalue problem, we restrict to functions ci 
and Si that are linear combinations of 6{p + p) and 6{p — p). This yields a four 
dimensional problem, whose eigenvalues can be calculated explicitly as functions 
of p, and therefore as functions of the energy per particle U . 

One finds that above U* = 7/12 all eigenvalues are purely imaginary, indicating that 
the Vlasov equation is linearly stable for such perturbations, in full agreement with 
the stability analysis developed in Sec. 3.4. Below U*, the largest eigenvalue A that 
controls the growth rate is 

A = ^QiU* - U) . (49) 

The exponential growth of perturbations of the initial distribution implies an expo- 
nential growth of magnetization M{t). This is indeed confirmed in Fig. 3(a) where 
the magnetization is plotted versus time in semi-logarithmic scale. The compari- 
son of the theoretical estimate (49) of the growth rate with the numerical results, 
reported in Fig. 3(b), shows a very good agreement. Moreover, the time scale 1/A 
diverges at f/*. 
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(b) 
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Fig. 3. (a) Semi-log plot of the magnetization M{t) for = 10^ particles, the number of 
samples being 20. The values of U are from 0.50 to 0.60 from top to bottom with 0.02 step 
size. The exponential growth rates of M{t), A, are estimated from the fitting solid lines, 
(b) Comparison of the theoretical (solid curve) and numerical (crosses) growth rate A. 



4-3 Extended initial distribution 



We present in this Section the numerical resuhs for the mixed distributions in for- 
mula (37) using the indicators introduced in Section 4.2.1. 

The first peak height, shown in Fig. 4, shows a dependence on U which is perfectly 
consistent with the theoretical prediction of the existence of a critical energy density 
U*{a), given in formula (38), above which the M = state is stable. 

However, the numerical results emphasize that the transition is much less abrupt 
for non vanishing values of the parameter a than for a = (the waterbag case 
previously analyzed). The explanation of this effect is possibly twofold. On one 
hand, the stability criterion derived in Sec. 3.3 gives no information on the time 
evolution that begins from an unstable initial state: it may well be that the system 
evolves initially to states with smaller magnetization than for a = 0. On the other 
hand, this weaker instability below U* may be also due to the linear vanishing of the 
growth rate at U* characteristic of the Gaussian initial distribution [4], instead of 
the sharper square root behaviour of the waterbag initial distribution, as expressed 
by Eq. (49). 

The first peak time data, presented in Fig. 5 as a function of the energy density U, 
do not follow sometimes a smooth curve. This may be due to the fact that magne- 
tization M{t) is almost flat in the region U ^ U*, and hence the first peak time is 
strongly affected by slight variations in the shape of the function M{t). Neverthe- 
less, the behavior of the first peak time is qualitatively the same as the one shown 
in Fig. 2(b) for a = 0. The same comments made there apply also to this case. 

We have shown in this Section that the formal stability criterion of Vlasov stationary 
states, stated in Section 3.3, allows to characterize the short time (0(1) time scale) 
behavior of the finite dynamics of the HMF model in the N oo limit. In 
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Fig. 4. The first peak heights versus the energy density U for waterbag (panel (a)), mixed 
with a = 0.5 (panel (b)) and Gaussian (panel (c)) initial conditions. The vertical line 
represents the theoretical critical energy density U*{a), given by Eq. (38). The different 
curves correspond to the following values of A^: 10^, 10'^, 10^ and 10^ from top to bottom 
for U > U*{a). 
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Fig. 5. The first peak time versus the energy density U for waterbag (panel (a)), mixed 
(panel (b)) and Gaussian (panel (c)) initial conditions. The vertical line represents the 
theoretical critical energy density U*{a), given by Eq. (38). The different curves correspond 
to the following values of N: 10^, 10'^, 10"^ and 10^ from bottom to top in U < U*{a). 



particular, it describes the behavior of the magnetization as a function of the energy 
per particle U and the presence of several different instability thresholds depending 
on the detailed properties of the chosen initial distribution. The critical energy for 
the instability of the Gaussian distribution coincides with the thermodynamical 
phase transition point Uc, as already remarked by Inagaki [4]. Moreover, the linear 
stability analysis of the Vlasov equation gives results for the growth rates of the 
instability that are perfectly confirmed by numerical simulations. Several aspects 
remain to be ascertained. An important one concerns the intermediate and long 
time evolution of an initial Vlasov unstable state. In the following Section, after 
introducing suitable indicators of stationarity we will follow the time evolution of 
such states for the HMF model. A further important question is the ultimate fate of 
Vlasov stable M = states for the finite dynamics of the HMF model. In the next 
Section, we will show that these states display indeed a much slower "instability", 
that becomes evident on a time scale that increases with a nontrivial power law 
in A^. These are examples of the quasi- stationary states that are the main object of 
study of this paper. 
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5 Intermediate and long time behaviours 

This Section is devoted to the study of the intermediate and long time evolution of 
both initially unstable and stable Vlasov stationary states. As discussed in the pre- 
vious section, the initial evolution is well described by Vlasov dynamics. Unstable 
states quickly evolve (on times of order 1), until they approach, and are trapped, 
close to stationary states of the Vlasov equation. The system then evolves slowly 
through these stationary states, until it reaches Boltzmann-Gibbs statistical equi- 
librium. In this section, we will construct indicators to assess if a particle state 
is close to a stationary state (necessary conditions). Using these indicators, we will 
show that in all cases, during the relaxation towards equilibrium, the system always 
remains close to Vlasov stationary states. 



5. 1 Necessary conditions of stationarity 

In order to check these features, we need to introduce necessary conditions for sta- 
tionarity, and to perform numerical tests of non stationarity for the HMF model. 
The tests guarantee non stationarity only, they do not guarantee stationarity. In 
order to obtain supporting evidences of stationarity, we introduce several indicators 
of non stationarity and observe whether all of them fail. We take this fact as a good 
indication of stationarity. 



5.1.1 Energy distribution 

If fst{diP) is stationary, so will be the single particle energy distribution, /e(e). This 
implies the stationarity of its moments 

^e^) = ^°°eVe(e)de , with j e N . (50) 

The logical implications 

f{6,p,t) stationary ^> feie,t) stationary ^> ^e-'^ (t) stationary (51) 

guarantees that the stationarity of the energy moments is a necessary condition for 
the stationarity of the distribution f{6,p,t). 

In other words, the system is not stationary if at least one of the moments (e-') (t) 
is not. The computation of the time derivatives d (e-^) /dt for the first four moments 
will be our Test I for non- stationarity. 
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5.1.2 Symmetry with respect to the spatial variable 6 



The previous necessary condition (51) is valid for all Hamiltonian models. It is 
possible to derive, using symmetry arguments, another necessary condition which 
is specific of the HMF model. Indeed, using the individual particle energy (14), the 
Vlasov equation (10) for a stationary distribution can be written as 



de d de d 
dp 09 86 dp 



fst{9,p) = . (52) 



The solution of Eq. (52) is constant on the characteristic curves e{6,p) =constant. 
Moreover, if one always resets to zero the phase of the magnetization vector M{t) 
(i.e. My{t) = 0), all characteristic curves will be symmetric with respect to a sign 
reversal of 6. Hence, if / is stationary and the previous condition on the phase of 
M{t) is respected, then the stationary distribution obeys the following symmetry 
condition 

fst{0,p) = fst{-0,p), \fp. (53) 
Consequently, an asymmetry of the distribution f{6,p,t) with respect to the spatial 
variable 6, implies that the system is non stationary. 

It is important to notice that the symmetry with respect to sign changes in mo- 
mentum p, f{9,p, t) = f{9, —p, t), is not a necessary consequence of the stationarity 
of /. Two separate characteristic curves for e > e<j, where Cg is energy of the sep- 
aratrix, may correspond to different values of / since they are not connected. This 
symmetry is therefore not required even if the system is stationary. A violation of 
this symmetry is not therefore an acceptable test of non stationarity. 

Obtaining the distribution f{6,p,t) numerically is, however, not easy since it is 
defined on the 2-dimensional /^-space at all fixed times t, and therefore a huge 
number of samples would be necessary to obtain a good statistics. A trick to easily 
check asymmetries in f{6,p,t) would be to verify the associated symmetry of the 
marginal distribution function 

f{9,t) = J f{9,p,t)dp . (54) 

However, unfortunately, f{0,t) may be symmetric even if the symmetry (53) is not 
satisfied. Typical examples of asymmetric distributions f{6,p,t) and corresponding 
symmetric marginal distributions f{6, t) are shown in Figs. 6(a) and (c), respectively. 
In order to recognize the asymmetry of f{0,p, t) with respect to 6, one must consider 
partially integrated distributions 

f+{9,t)=f f{9,p,t)dp and U9,t) = f f{e,p,t)dp . (55) 

Jp>0 Jp<0 

These latter distributions, besides being symmetric if the system is stationary, are 
able to correctly detect the asymmetries of 2D distributions in {6,p), as shown in 
Figs. 6(a) and (c). 
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Fig. 6. Typical asymmetric (a) and symmetric (b) distribution f(9,p,t). Panel (c) 
(resp. (d)) presents the marginal distributions f{9,t) (crosses), f+{6,t) (open squares) 
and f-{6,t) (full squares) for the distributions plotted in panel (a) (resp. (b)). The three 
curves superpose almost perfectly in panel (d). = 10^, U = 0.63, a = 1.0, the number 
of samples is 10^. Time is t = 10 for panels (a) and (c), and t = 10^ for panels (b) and 
(d). 
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We quantitatively observe the time evolution of the asymmetry by introducing the 
two following quantities 



A4t) 



U{9,t) - U{-9,t) 

Uo,t) - f4-e,t) 



de 



d9 



de 
de 



p>0 



p<0 



dp{f{e,p,t)-f{-e,p,t)) 
dpifie,p,t)-fi-e,p,t)) 



(56) 



Both y4_|_ and A_ exactly vanish when the distribution is symmetric, but they are 
affected by numerical errors and finite effects when estimated from numerical 
simulations. In order to reduce such effects, we first define the quantity 



A{t) 



de 



dp{f{e,p,t)-f{-e, -p,t)) 



(57) 



For initial conditions that are symmetric under the transformation (6*, p) (— ^, —p), 
since this symmetry is conserved during the time evolution, A(t) = for both sta- 
tionary or non stationary distributions. Of course, also this quantity is affected by 
numerical errors. To factor out such errors, we introduce the quantities A+{t)/A{t) 
and A_{t) /A(t), which are large for asymmetric distributions similar to those in 
Fig. 6(a), and are instead 0(1) for symmetric distributions (actual values taken for 
the cases shown Fig. 6 are 160 for Fig. 6(a) and 0.4 for Fig. 6(b)). Large values of 
the quantities A±/A will be our Test II for non stationarity. 



5.2 Non- stationarity test 



Now that we have defined these two indicators, let us check non stationarity by 
considering both Vlasov unstable and Vlasov stable initial states with waterbag 
distribution, a = in Eq.(37). Typical results of Test-I and -II are shown in Fig. 7. 

For Vlasov unstable cases (Fig. 7(a)), the quantity d (e-') /dt takes large values only 
in the time region 2 < t < 80. Similarly, the indicator A±{t)/A{t) is large in the 
time region 1 < t < 10. We have reported only d (e) /dt and Aj^{t)/A{t), but we 
have checked that the behaviour is totally identical for higher moments j = 2, 3 and 
4 and for A_{t) / A{t), respectively. This time region of non stationarity perfectly 
coincides with the region where the magnetization M{t) rapidly grows or largely 
fluctuates. For Vlasov stable cases, d (e-') /dt exhibits small fluctuations only in the 
time interval 1 < t < 20, while A±{t)/A{t) remains small all the time. We guess that 
the fluctuations in d(e-') /dt are due to finite size effects and would vanish in the 
N ^ oo limit. Note the logarithmic scale for the horizontal time axis: the increase 
in M{t) in Fig. 7(b) is much slower than the one in Fig. 7(a). 

The behavior shown in Fig. 7(b) is very important, because it shows that a variation 
in M is compatible with the fact that the distribution / may remain stationary; i.e. 
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Fig. 7. Non stationarity tests for a waterbag initial condition with N = 10^ particles 
and 10^ samples. Note the logarithmic scale for the time variable. Panels (a) (resp. (b)) 
presents the time evolution of the magnetization M{t) (crosses), Test I d (e) /dt (open 
squares) and Test II A^{t)/A{t) (full squares) in the unstable case U = 0.55 (resp. stable 
case U = 0.69). The quantities d^e-') /dt (resp. A±{t)/A{t)) are multiplied by a factor 
10 (resp. 10~^) for graphical purposes. The horizontal line represents the canonical value 
of M. 



the system can relax to equilibrium (notably in the time region [10^, 10^] in Fig. 7(b)) 
passing through stationary states of the Vlasov equation. It is natural to suppose in 
addition that these stationary states are stable. 

Figure 8 reports results for two unstable Gaussian initial conditions, a = 1 in 
Eq.(37)with U = 0.55 and U = 0.69. The behavior in panel (a) is very similar 
to the unstable waterbag case of Figs. 7(a). The case in panel (b) exhibits a small 
peak around t = 5 for both Test I and II. This presumably indicates that the sys- 
tem is weakly non stationary at this time, meaning that the temporal variations of 
the distributions / are small. It is also possible to explain the weakness of the non 
stationarity for Gaussian initial conditions with the same arguments that we used 
when commenting Fig. 4. This state is unstable below the critical energy Uc = 0.75, 
but the equilibrium value for M goes to zero when U tends to Uc- Hence, the initially 
unstable state may be very close to stationary stable states, and can be therefore 
immediately trapped into a stable state. We expect the non stationarity to become 
weaker and weaker as U approaches Uc- 

The initial conditions used above correspond to stationary solutions of the Vlasov 
equation. What happens with non stationary initial conditions ? On the basis of 
what we have found above, we expect that the system reaches a stationary stable 
state, just like it does when starting from initially unstable stationary states (Figs. 7 
and 8). To confirm this expectation, we have prepared the initial condition M(0) = 1, 
and we have reported the results of Tests I and II in Fig. 9. In the short time region 
t < 3, both d (e-') /dt and A± are large but they vanish for longer times, indicating 
that the system reaches some stationary state of the Vlasov equation. 

Summarizing the results of this Section, we have found a strong support to the 
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Fig. 8. Same plots as in Figs. 7, for a Gaussian initial condition. 
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Fig. 9. Same plots as in Figs. 7, for waterbag initial conditions in momentum p whereas 
all the phases 6j are set at 0, i.e. M(0) = 1. The quantity Aj^{t)/A{t) is multiplied by 
10~^ for graphical purposes. 

guess that, generically, the states corresponding to stationary stable solutions of 
the Vlasov equation attract stationary unstable and non stationary solutions. The 
system evolves in time passing through a series of stationary solutions of the Vlasov 
equation. 



5.3 Long time behaviour 



Let us now consider the long time behaviour of the system. We first discuss the be- 
havior of the time evolution of the marginal distribution in angle /(^, t) = J f{6, p, t)dp, 
shown in Fig. 10. For U = 0.55, the initial condition is Vlasov unstable, thus / dras- 
tically changes from t = 1 to t = 10 (see Figs. 10(a) and 7). From t = 10 to t = 10^, 
the distribution is almost frozen and finally relaxes to the equilibrium distribution 
at t = 10^. On the contrary, in the Vlasov stable case, U = 0.69, (see Fig. 10(b)) the 
distribution is almost frozen from t = 1 up to t = 10'^: this is the crucial test of the 
presence of a quasi- stationary state associated to the Vlasov stable initial condition. 
Meanwhile, the magnetization M{t) remains constant, well below the equilibrium 
value. Finally the distribution changes from t = 10^ to 10^, and correspondingly 
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the magnetization M{t) increases to reach equihbrium (similar results have been 
obtained for Gaussian initial distributions). 




Fig. 10. Temporal evolution of the probability distribution function f{9, t) for a waterbag 
initial condition, identical to the one used for Fig. 7. N = 10^ and the number of samples 
is 10^. The distributions at times t = I (crosses), t = 10 (open squares), t = 10^ (full 
squares), t = 10^ (open circles), are reported for U = 0.55 (a), and at t = 1 (crosses), 
t = 10^ (open squares), t = 10^ (full squares), t = 10^ (open circles), for U = 0.69 (b). The 
solid curve represents the equilibrium distribution, exp[M cos 0/T]/C, where M = 0.558, 
T = 0.412 and C = 9.52 for U = 0.55 and M = 0.309, T = 0.475 and C = 6.97 for 
U = 0.69. 

In order to characterize the relaxation timescale with respect to the particle num- 
ber in the case of the Vlasov stable waterbag initial conditions {U = 0.69), we 
have studied the temporal evolution of the magnetization M{t) (a preliminary study 
has been already reported in Ref. [34]). This can be fitted by the function 

M{t) = [1 +tanh(a(A^)(logiot- W))] c{N)+d{N), (58) 

as shown in Fig. 11(a). The parameters c{N) and d{N) represent the half width be- 
tween the initial and the equilibrium levels of M{t) and the initial level of M{t), re- 
spectively (we further comment about such parameters in the Appendix). The prod- 
uct a(A^)c(A^) is the slope at logio t = b{N),i.e. a{N)c{N) = dM/d(logiot)|iogioi=fe{A^)> 
and t{N) = 10^*^^-' is the timescale. The most important parameter, t{N), presented 
in Fig. 11(b) as a function of TV, is shown to be proportional to A^^''^. The fit is very 
good and excludes both the and the A^^ trivial scalings. This nontrivial power 
law emerges unexpectedly, since theoretical arguments as well as previous studies 
of the HMF model suggest trivial divergences of the timescale (typically as A^) (see 
Ref. [3]). Quite interestingly however, Zanette and Montemurro [25] have also in- 
vestigated the timescale for the HMF model for U = 0.69 with a waterbag initial 
distribution of momenta, but with a fully ordered initial state, i.e. M(0) = 1. They 
have shown that M{t) presents a minimum at tmin (already present in the simula- 
tion by Latora et al [11]) which is also proportional to A^^'^. Although our initial 
condition is different and hence different scalings could be found, we consider this 
finding as a confirmation of the presence of nontrivial timescales in the HMF model. 
Other, apparently unrelated, nontrivial timescales have been found, associated to 
the vanishing Lyapunov exponent, for instance the well verified 1/3 scaling [3]. 
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Fig. 11. Panel (a) presents the temporal evolution of the magnetization M(i) for different 
particles numbers: N = 102(10^), 10^(102), 2.10^(8), 5.10^(8), 10^(8) and 2.10-^(4) from left 
to right, the number between brackets corresponding to the number of samples. The hori- 
zontal line represents the equilibrium value of M. Panel (b) shows the logarithmic timescale 
b{N) as a function of N, whereas the dashed line represents the law lO''^^-* ~ A^-^''''. 

We were led by this discovery to rescale the momentum distributions f*{p,t) = 
I fi^^yPy i)dO by the nontrivial time-scale r(iV) ~ iV^-^, see Fig. 12. The superposition 
of the curves at different times, in the quasi-stationary state, is quite impressive and 
we consider this as a strong test of the nontrivial timescale. 




Fig. 12. Temporal evolution of the momentum distribution f*{p,t/T{N)) at times rescaled 
by iV^-'^ for the Vlasov stable case U = 0.69: N = 10^ (pluses), = 2.10^ (crosses), 
N = 3.10^ (stars) and N = 10^ (open squares). The distributions are plotted for 
t/N^'''' = 0.0095, for instance, t = 6.10^ for = 10^. The four curves, once time-scaled, 
superpose almost perfectly. 

There have been recently claims that momentum distributions in quasi-stationary 
states, obtained from the M(0) = 1 initial conditions, could be fitted with single 
particle distributions inspired by Tsallis statistics [11]. In this study, we could not 
obtain any reasonable fit of the momentum distribution in Fig. 12 using several 
values of the Tsallis q parameter with g > 1. Moreover, the tails of the distribution 
are rather sharp, as shown by the inset in lin-log scale, and we can numerically 
exclude the presence of power law tails. 



It should be however mentioned that in another physical context, i.e. in the dy- 
namics of self-gravitating A^-body systems confined in an adiabatic wall, the long 
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term time evolution has been argued to sweep through stellar polytropes, that are 
pecuhar stationary solutions that maximize Tsallis entropy [35]. However, Chava- 
nis [36], using the minimization of functionals built on conserved quantities (called 
H-functions in this context), showed that Tsallis entropy can be considered as a par- 
ticular choice among the infinitely many possible H-functions: their maximization 
lead to particular stationary solutions of the Vlasov dynamics. 

The success of Tsallis statistics in describing quasi- stationary regimes in N-body 
dynamics remains therefore doubtful. Further studies should fix this puzzling issue. 

It should be moreover mentioned that in connection with some aging properties of 
the M = 1 initial state, g-exponential have been used to fit quite efficiently numerical 
data of correlation functions [37,13]. 



6 Summary 

In this paper, we have proposed a general framework, based on Vlasov equation, to 
study the relaxation to equilibrium of the HMF model. In the short time regime, the 
behaviour of the system is well described by the Vlasov equation, which is a valid 
approximation for finite times. The stability of the initial conditions is therefore 
adequately characterized by considering the stability of Vlasov stationary states. We 
have checked the accuracy of this criterion in comparison with finite simulations. 
In addition, a simple analytical derivation of the largest eigenvalue of the linearized 
Vlasov operator gives an excellent estimate of the numerical timescale in the short 
time regime of the finite system. 

In the intermediate and long time regimes, two situations have been analyzed sep- 
arately. The first one, corresponding to unstable stationary Vlasov states, shows a 
complex short time evolution before reaching some stationary stable state of the 
Vlasov equation that is situated close to Boltzmann-Gibbs equilibrium. The latter 
is finally reached on very long times, t{N) ~ A^. The second one, which begins from 
Vlasov stable stationary states, evolves through other quasi- stationary states that 
are far from equilibrium and are among the many other stationary Vlasov states. The 
final relaxation to Boltzmann-Gibbs equilibrium takes place on times that increase 
with a nontrivial power law in the number of particles t{N) ~ A^^'^. Interestingly, 
this exponent is consistent with the one recently reported in Ref. [25] for a different 
initial condition. 

We have analyzed the momentum distributions in such a quasi- stationary regime 
and found that they scale properly with the nontrivial exponent. As in [13], we 
have been unable to fit the momentum distributions using g-exponentials, as it had 
instead been done for other initial states in Ref. [11]. This is a negative indication 
that quasi- stationary states can always be described by Tsallis statistics. 
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An analogy between the slow time evolution appearing in quasi- stationary states and 
the aging phenomenon in glasses and spin-glasses has been proposed in Ref. [37] and 
further analyzed in Ref. [13]. However, none of these studies has been performed 
using stationary states of the Vlasov equation as initial conditions. It would be 
extremely interesting to repeat the study of correlation functions and of other glassy 
behavior indicators for such states, as already partially done in Ref. [34] , where non 
stationary stretched exponential correlation functions have been found. 

After completing this paper, we became aware of a further, and more complete, 
study by Choi and Choi [38] of the linear stability of the homogeneous state in both 
the canonical (Fokker- Planck) and the microcanonical (Vlasov) ensemble. However, 
this paper does not give any answer to the question of the instabilities caused by 
the finiteness of the number of particles. 
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A Appendix: Detailed scaling of the magnetization 



Let us recall the fitting function we have used in Section 5.3 for the magnetiza- 
tion M{t) 



We have already commented about the scaling behavior of t{N) = lO^*-^-*. According 
to Fig. A.l, all the parameter set scales as follows: 



M{t) = [l + tanh(a(iV)(logiot-&(iV)))] c{N) + d{N) . 



(A.l) 
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where M^q is the equilibrium level of M{t). 
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Fig. A.l. Log-log plots of parameters a{N)c{N) and d{N) in panel (a), and c{N) in panel 
(b). The symbols represent numerical results whereas the curves correspond to the fitted 
functions, given by Eqs. (A. 2). 

By using the scaling functions (A. 2), we can predict when M{t) reaches a given 
threshold, Mth, as a function of A^. The corresponding threshold time, tth, with 
^(^th) = Mth, has the following expression 
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Let us consider, for numerical purposes, two threshold times for both threshold levels, 
Mth = d+e and Meq— These two times roughly represent the beginning and ending 
times when M{t) grows towards the equilibrium value Meq [34]. Numerical results 
have fluctuations, particularly in the early time region, and hence we have defined 
tth{d + e) and tth (Meq — as follows: 



tth{d + e)=max{t e R\ M{t) <d + e} , (A.4) 
tth(Meq -e)= mm{t e R \ M{t) > Meq - e} . (A.5) 

In Fig. A. 2, the comparison between the numerical results and the theoretical pre- 
diction (A.3) for e = 0.0088 shows an excellent agreement, except in the small A^ 
region. We remark that both threshold times are approximately proportional to A^^ '' 
in the thermodynamic limit. This further confirms the nontrivial scaling law. 
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